#include "VGQApprox.h"
#include <fstream>
#include <algorithm>
#include <boost\math\tools\test_data.hpp>
#include "CFraction.h"


VGQApprox::VGQApprox( double mu, double lambda, double alpha, double beta, double eps )
	: _mu(mu), _pvgDistribution(new VarianceGammaDistribution(0.0, lambda, alpha, beta)), 
		_asymptoticExpansion(mu, lambda, alpha, beta), _eps(eps), _maxDegree(50), _maxAccuracy(0.0)
{
	initialize();
}

VGQApprox::~VGQApprox(void)
{
}

void VGQApprox::initialize()
{
	using namespace std;

	_pcentralLSeries = make_shared<VGQSeries>(_pvgDistribution);
	_pcentralRSeries = make_shared<VGQSeries>(_pvgDistribution);
	_pleftSeries = make_shared<VGQLeftSeries>(_pvgDistribution);
	_prightSeries = make_shared<VGQRightSeries>(_pvgDistribution);

	//VGQSeries testSeries(_pvgDistribution);
	//testSeries.SetInitialCondition(0.5);
	//testSeries.radius();

	// determine the cut off points
	_u0 = _pvgDistribution->zeroQuantileLocation();
	findLeftTailCutOff();
	findRightTailCutOff();

	// Build Approximants
	buildCentralLeftApprox();
	buildCentralRightApprox();
	buildLeftApprox();
	buildRightApprox();

}

void VGQApprox::buildCentralRightApprox()
{
	double u0CR = (1 + 4*_u0)/5;

	for(int k = 0; k <= 50; k++)
	{
		_pcentralRSeries->SetInitialCondition(u0CR);
		_uR = min(u0CR + _pcentralRSeries->radius(), 1 - std::numeric_limits<double>::epsilon());
		
		double x0CR = _pcentralRSeries->GetInitialCondition().second;
		_approxCR = make_shared<CFraction>(_pcentralRSeries, x0CR);

		// Build central left approximation
		for(unsigned int n = 6; n <= _maxDegree; n += 2)
		{
			_approxCR->build(n);
			_approxCR->convert();

			double q1 = x0CR + _approxCR->evaluate(_u0 - u0CR);
			double q2 = x0CR + _approxCR->evaluate(_uR - u0CR);

			if(uError(_u0, q1) < _eps && uError(_uR, q2) < _eps)
				return;
		}

		u0CR = (_u0 + u0CR)/2;

	}

	throw std::runtime_error("Unable to satisfy accuracy requirements in right central region.");
}

void VGQApprox::buildCentralLeftApprox()
{
	double u0CL = (4*_u0)/5;
	
	for(int k = 1; k <= 50; k++)
	{
		_pcentralLSeries->SetInitialCondition(u0CL);
		//_uL = max(_u0 - _pcentralLSeries->radius(), 0.1);
		_uL = max(u0CL - _pcentralLSeries->radius(), std::numeric_limits<double>::epsilon());
		double x0CL = _pcentralLSeries->GetInitialCondition().second;
		_approxCL = make_shared<CFraction>(_pcentralLSeries, x0CL);

		// Build central left approximation
		for(unsigned int n = 6; n <= _maxDegree; n += 2)
		{
			_approxCL->build(n);
			_approxCL->convert();

			double q1 = x0CL + _approxCL->evaluate(_uL - u0CL);
			double q2 = x0CL + _approxCL->evaluate(_u0 - u0CL);

			if(uError(_uL, q1) < _eps && uError(_u0, q2) < _eps)
				return;
		}

		u0CL = (_u0 + u0CL)/2;
	}
	
	throw std::runtime_error("Unable to satisfy accuracy requirements in left central region.");

}

void VGQApprox::buildRightApprox()
{
	double step = 0.1;
	double t = 0.9;

	// for t = 0.9, 0.8, ..., 0.1
	for(int k = 0; k <= 8; k++)
	{
		t -= step;

		// Set the initial conditions
		double u0R = ((1-t) * _uR) + (t * _tauR);
		_prightSeries->SetInitialCondition(u0R);

		double x0R = _prightSeries->GetInitialCondition().second;
		double z0R = _prightSeries->GetInitialCondition().first;
		_approxR = make_shared<CFraction>(_prightSeries, x0R);

		// Build left approximation 
		for(unsigned int n = 6; n <= _maxDegree; n += 2)
		{
			_approxR->build(n);
			_approxR->convert();

			double q1 = x0R + _approxR->evaluate(_prightSeries->baseQuantile(_uR)-z0R);
			double q2 = x0R + _approxR->evaluate(_prightSeries->baseQuantile(_tauR)-z0R);

			if(uError(_uR, q1) < _eps && uError(_tauR, q2) < _eps)
				return;
		}
	}

	throw std::runtime_error("Unable to satisfy accuracy requirements in right region.");
}

void VGQApprox::buildLeftApprox()
{
	double step = 0.1;
	double t = 0.9;

	// for t = 0.9, 0.8, ..., 0.1
	for(int k = 0; k <= 8; k++)
	{
		t -= step;

		// Set the initial conditions
		double u0L = (t * _tauL) + ((1-t) * _uL);
		_pleftSeries->SetInitialCondition(u0L);
		
		double x0L = _pleftSeries->GetInitialCondition().second;
		double z0L = _pleftSeries->GetInitialCondition().first;
		_approxL = make_shared<CFraction>(_pleftSeries, x0L);

		// Build left approximation 
		for(unsigned int n = 6; n <= _maxDegree; n += 2)
		{
			_approxL->build(n);
			_approxL->convert();

			double q1 = x0L + _approxL->evaluate(_pleftSeries->baseQuantile(_uL)-z0L);
			double q2 = x0L + _approxL->evaluate(_pleftSeries->baseQuantile(_tauL)-z0L);

			if(uError(_uL, q1) < _eps && uError(_tauL, q2) < _eps){	
				return;
			}
		}
	}

	throw std::runtime_error("Unable to satisfy accuracy requirements in left region.");
}

double VGQApprox::operator()( double u )
{
	if(0 <= u && u < _tauL)
		return _asymptoticExpansion.LeftTail(u);

	if(_tauL <= u && u < _uL){
		double z0 = _pleftSeries->GetInitialCondition().first;
		double x0 =_pleftSeries->GetInitialCondition().second;
		return _mu + x0 + _approxL->evaluate(_pleftSeries->baseQuantile(u)-z0);
	}

	if(_uL <= u && u < _u0){
		double p0 = _pcentralLSeries->GetInitialCondition().first;
		double x0 =_pcentralLSeries->GetInitialCondition().second;
		return _mu + x0 + _approxCL->evaluate(u-p0);
	}

	if(_u0 <= u && u < _uR){
		double p0 = _pcentralRSeries->GetInitialCondition().first;
		double x0 =_pcentralRSeries->GetInitialCondition().second;
		return _mu + x0 + _approxCR->evaluate(u-p0);
	}

	if(_uR <= u && u <= _tauR){
		double z0 = _prightSeries->GetInitialCondition().first;
		double x0 =_prightSeries->GetInitialCondition().second;
		return _mu + x0 + _approxR->evaluate(_prightSeries->baseQuantile(u)-z0);
	}

	if(_tauR < u && u <= 1)
		return _asymptoticExpansion.RightTail(u);

	throw std::domain_error("u not in [0,1]");
}

double VGQApprox::uError( double u, double xApprox )
{
	double uRes = std::abs(u - _pvgDistribution->cdf(xApprox));
	_maxAccuracy = std::max(uRes, _maxAccuracy);
	return uRes;
}

void VGQApprox::profile( const std::string& directoryPath )
{
	using namespace boost::math::tools;
	int nPoints = 25;
	double stepSize = _tauL / nPoints;
	
	// left tail region
	std::string fileName = directoryPath + "\\leftTail.csv";
	shared_ptr<ofstream> fleftTail = openFile(fileName);

	test_data<double> testData;
	testData.insert(*this, make_periodic_param(0 + stepSize, _tauL, nPoints));
	write_csv(*fleftTail, testData, ", ");

	// left region
	fileName = directoryPath + "\\left.csv";
	shared_ptr<ofstream> fleft = openFile(fileName);

	nPoints = 100;
	stepSize = (_uL - _tauL) / nPoints;
	testData.clear();
	testData.insert(*this, make_periodic_param(_tauL + stepSize, _uL, nPoints));
	write_csv(*fleft, testData, ", ");

	// left central region
	fileName = directoryPath + "\\leftCentral.csv";
	shared_ptr<ofstream> fleftCentral = openFile(fileName);

	nPoints = 100;
	stepSize = (_u0 - _uL) / nPoints;
	testData.clear();
	testData.insert(*this, make_periodic_param(_uL + stepSize, _u0, nPoints));
	write_csv(*fleftCentral, testData, ", ");

	// right central region
	fileName = directoryPath + "\\rightCentral.csv";
	shared_ptr<ofstream> frightCentral = openFile(fileName);

	nPoints = 100;
	stepSize = (_uR - _u0) / nPoints;
	testData.clear();
	testData.insert(*this, make_periodic_param(_u0 + stepSize, _uR, nPoints));
	write_csv(*frightCentral, testData, ", ");

	// right region
	fileName = directoryPath + "\\right.csv";
	shared_ptr<ofstream> fright = openFile(fileName);

	nPoints = 100;
	stepSize = (_tauR - _uR) / nPoints;
	testData.clear();
	testData.insert(*this, make_periodic_param(_uR + stepSize, _tauR, nPoints));
	write_csv(*fright, testData, ", ");

	// right tail region
	fileName = directoryPath + "\\rightTail.csv";
	shared_ptr<ofstream> frightTail = openFile(fileName);
	
	// Special handling for the right tail (for some reason testData.insert does not work in this case)
	nPoints = 25;
	stepSize = (1 - _tauR) / nPoints;
	double a[25];
	for(int n = 0; n < 24; n++)
	{
		a[n]=_tauR + stepSize*double(n+1);
	}
	for(int n = 0; n < 24; n++)
	{
		*frightTail << a[n] << ", " << (*this)(a[n]) << endl;
	}
	
	//testData.clear();
	//testData.insert(*this,make_periodic_param(_tauR, 1.0, nPoints));
	//write_csv(*frightTail, testData, ", ");

	cout << "Left Approximant Degree: " << _approxL->degree() << endl;
	cout << "Left Central Approximant Degree: " << _approxCL->degree() << endl;
	cout << "Right Central Approximant Degree: " << _approxCR->degree() << endl;
	cout << "Right Approximant Degree: " << _approxR->degree() << endl;
}

shared_ptr<ofstream> VGQApprox::openFile( std::string fileName )
{
	shared_ptr<ofstream> pFile = make_shared<ofstream>(fileName);
	if(!pFile->good())
		throw runtime_error(std::string("Failed to open file: ") + fileName);

	pFile->precision(16);
	return pFile;
}

void VGQApprox::findLeftTailCutOff()
{
	// We want to find a region (0, _taulL] for which the asymptotic expansion is valid
	double xL = _pvgDistribution->getSupport().first;
	double minCutOff = std::numeric_limits<double>::epsilon();
	_tauL = std::max(minCutOff, _pvgDistribution->cdf(xL));

	int k = 0;
	do
	{		
		if(uError(_tauL, _asymptoticExpansion.LeftTail(_tauL) - _mu) < _eps)
			return;

		if(_tauL <= minCutOff)
			break;

		k++;
		_tauL /= 2;
	} while(k < 50);

	throw std::runtime_error("Unable to satisfy accuracy requirements in left tail region.");
}

void VGQApprox::findRightTailCutOff()
{
	// We want to find a region [_taulR,1) for which the asymptotic expansion is valid
	double xR = _pvgDistribution->getSupport().second;
	double maxCutoff = 1 - std::numeric_limits<double>::epsilon();
	_tauR = std::min(maxCutoff, _pvgDistribution->cdf(xR));

	int k = 0;
	do
	{		
		if(uError(_tauR, _asymptoticExpansion.RightTail(_tauR) - _mu) < _eps)
			return;

		if(_tauR >= maxCutoff)
			break;

		k++;
		_tauR = std::min(maxCutoff, (1+_tauR)/2);
	} while(k<50);

	throw std::runtime_error("Unable to satisfy accuracy requirements in right tail region.");
}




